Code
library(tidyverse)
library(sf)
library(tmap)
library(terra)
library(stars)
library(kableExtra)December 3, 2025
Link to GitHub repository: https://github.com/aakriti-poudel-chhetri/impacts-of-extreme-weather
From February 13 to 17, 2021, Winter Storm Uri brought unusually cold temperatures, heavy snowfall, ice and freezing rain across Texas. Many places experienced their coldest conditions in more than three decades. The state’s power infrastructure was not designed for such extreme cold, therefore numerous power plants and grid components froze or malfunctioned. At the same time, electricity demand surged as residents attempted to heat their homes, placing severe pressure on the power system. To avoid a complete grid failure, authorities implemented rolling blackouts that were intended to last about an hour, but in reality many households lost electricity for several hours or even longer. An estimated 4.5 million homes experienced outages. These disruptions triggered additional problems, including widespread water and food shortages, burst pipes, and a lack of indoor heating. At least 246 deaths were linked directly or indirectly to the storm. Economic losses exceeded $195 billion, making Winter Storm Uri the expensive natural disaster in Texas history(Zhou et al. 2024).
Understanding how communities were affected is crucial for preparing for future extreme weather events. It can help identify which neighborhoods are most vulnerable and guide decisions about where to focus recovery and emergency resources.
Thus, this study examines how Winter Storm Uri affected communities across the Houston metropolitan area by estimating the number of homes that lost power and assessing whether lower-income neighborhoods faced a disproportionate share of the impacts. Using satellite-based nighttime lights data from the VIIRS VNP46A1 product, we identify areas where electricity was likely lost by comparing light intensity before and after the storm. These blackout areas are then combined with OpenStreetMap building data to estimate the number of affected homes and linked with U.S. Census Bureau data to evaluate socioeconomic differences—such as variations in median household income—between impacted and non-impacted census tracts.
Understanding how power outages unfolded across Houston during Winter Storm Uri is important for evaluating both the extent of the blackout and the communities most heavily affected. Using satellite-based detection of outage areas alongside socioeconomic data, this study examines whether the storm’s impacts were distributed unevenly across neighborhoods. The main research question guiding this analysis is:
“How were power outages during Winter Storm Uri distributed across the Houston metropolitan area, and did these outages disproportionately affect lower-income neighborhoods?”
This study also addresses the following sub-questions:
For this analysis, data were extracted from the multiple sources.
VIIRS data is distributed through NASA’s Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC). Many NASA Earth data products are distributed in 10x10 degree tiles in sinusoidal equal-area projection. Tiles are identified by their horizontal and vertical position in the grid. Houston lies on the border of tiles h08v05 and h08v06. These two tiles per date were pre-downloaded and provided by the team.
To assess the impacts of Winter Storm Uri on Houston communities, a multi-step approach combining spatial and socioeconomic data was adopted. This allowed us to identify areas affected by power outages, estimate the number of homes impacted, and examine patterns across different socioeconomic groups. The major steps of the analysis are summarized below:
Let’s start by importing all the necessary libraries.
First, we read our datasets.
# Read night lights data
tile05_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v05.001.2021039064328.tif'))
tile06_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v06.001.2021039064329.tif'))
tile05_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather','data', 'VNP46A1', 'VNP46A1.A2021047.h08v05.001.2021048091106.tif'))
tile06_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021047.h08v06.001.2021048091105.tif'))
# Read roads data
highways <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'gis_osm_roads_free_1.gpkg'),
query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass = 'motorway'") %>%
st_transform(crs = 'EPSG:3083')
# Read houses data
houses <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data',
'gis_osm_buildings_a_free_1.gpkg'),
query = "SELECT *
FROM gis_osm_buildings_a_free_1
WHERE (type IS NULL AND name IS NULL)
OR type IN ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>%
st_transform('EPSG:3083')
# Explore the contents of the geodatabase
socioeconomic <- st_layers(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'))
# Extract the census tract information
census_tract <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'ACS_2019_5YR_TRACT_48_TEXAS') %>%
st_transform('EPSG:3083')
# Extract the income information
income <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'X19_INCOME')We will merge the raster tiles and crop it to the Houston bounding box to focus the analysis on the study area. It will ensure only valid night-light data within Houston were used for assessing power outages.
# Merge two raster data of February 7
nightlight_07 <- st_mosaic(tile05_07, tile06_07)
# Merge two raster data of February 16
nightlight_16 <- st_mosaic(tile05_16, tile06_16)
# Create Houston bounding box
houston_bbox <- st_bbox(c(xmin = -96.5, xmax = -94.5, ymin = 29, ymax = 30.5),
crs = 'EPSG:4326')
# Crop all raster data to Houston area
nightlight_07_crop <- st_crop(nightlight_07, houston_bbox)
nightlight_07_crop[(nightlight_07_crop > 1000) | (nightlight_07_crop <= 0)] <- NA
nightlight_16_crop <- st_crop(nightlight_16, houston_bbox)
nightlight_16_crop[(nightlight_16_crop > 1000) | (nightlight_16_crop <= 0)] <- NAHere, we visualize the comparison between night-light intensities in Houston before and after the February 2021 winter storm. By mapping the February 7 (pre-storm) and February 16 (post-storm) night-light data side by side, the analysis shows areas where lighting and likely power availability changed, helping to identify regions affected by the blackout.
tmap_options(component.autoscale = FALSE)
# Map for February 7 (before storm)
beforestorm <- tm_shape(nightlight_07_crop) +
tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
col.legend = tm_legend(show = TRUE,
title = 'Night light intensity before storm',
title.size = 1,
orientation = 'landscape',
width = 35,
height = 3,
legend.just = c(0, 1)))+
tm_title('Houston Night Lights\nFebruary 07, 2021 (Before storm)',
size = 1.2, fontface = 'bold') +
tm_compass(position = c('left', 'bottom'),
type = 'arrow',
size = 1.5,
text.size = 0.7,
color.dark = 'grey50',
text.color = 'white') +
tm_scalebar(position = c('left', 'bottom'),
text.size = 0.8,
color.dark = 'grey50',
text.color = 'white') +
tm_layout(bg.color = 'white',
outer.bg.color = 'white',
frame = TRUE,
legend.show = TRUE) +
tm_graticules(x.show = TRUE, y.show = FALSE, lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)
#tm_graticules(lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)
# Map for February 16 (after storm)
afterstorm <- tm_shape(nightlight_16_crop) +
tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
col.legend = tm_legend(show = TRUE,
title = 'Night light intensity after storm',
title.size = 1,
orientation = 'landscape',
width = 35,
height = 3)) +
tm_title('Houston Night Lights\nFebruary 16, 2021 (After storm)',
size = 1.2, fontface = 'bold') +
tm_compass(position = c('left', 'bottom'),
type = 'arrow',
size = 2.5,
text.size = 0.7,
color.dark = 'grey50',
text.color = 'white') +
tm_scalebar(position = c('left', 'bottom'),
text.size = 0.8,
color.dark = 'grey50',
text.color = 'white') +
tm_layout(bg.color = 'white',
outer.bg.color = 'white',
frame = TRUE,
legend.show = TRUE) +
tm_graticules(x.show = TRUE, y.show = FALSE, lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)
tmap_arrange(beforestorm, afterstorm, ncol = 2)We will now create a blackout mask by identifying areas of Houston where night light intensity decreased significantly after the storm, indicating potential power outages. By calculating the difference between before and after night light, cropping to the study area, filtering out minor changes, and vectorizing the result, the mask provides a spatial representation of blackout affected regions for further analysis.
# Calculate the difference in light intensity
nightlight_diff <- nightlight_07 - nightlight_16
# Crop and reclassify the raster data
mask <- st_crop(nightlight_diff, houston_bbox)
mask[mask < 200] <- NA
# Vectorize the blackout mask
blackout <- st_as_sf(mask,
as_points = FALSE,
merge = TRUE) %>%
st_make_valid() %>%
st_transform(crs = 'EPSG:3083')In this step, highway geometries are combined and buffered to exclude areas where changes in night light intensity could be due to traffic rather than power outages. By removing blackout areas within 200 meters of highways, the analysis isolates regions where lighting loss more likely reflects true residential or commercial power outages.
# Combine all highway geometries into one
highways_union <- st_union(highways)
# Create 200m buffer around all highways
highways_buffer <- st_buffer(highways_union, dist = 200)
# Find areas that experienced blackouts further than 200m from highways
blackout_far <- st_difference(blackout, highways_buffer)We visualize the spatial extent of power outages across Houston by mapping the blackout areas. The map highlights the neighborhoods affected, providing a clear view of the geographic distribution of the electricity disruption during the storm.
# Visualize the blackout in Houston
tm_basemap('OpenStreetMap') +
tm_shape(blackout_far) +
tm_polygons(fill = 'ivory', fill_alpha = 0.8, col = 'black') +
tm_title('Blackout areas in Houston', fontface = 'bold', size = 1) +
tm_layout(legend.show = TRUE) +
tm_scalebar(position = c('left', 'bottom')) +
tm_compass(position = c('right', 'top')) +
tm_graticules(lwd = 0.2, alpha = 0.3, labels.size = 0.5) +
tm_add_legend(labels = c("Blackout Areas"),
fill = c('black'),
type = 'polygons',
position = c('left', 'top'))We estimate the number of homes in Houston that lost power during the February 2021 winter storms. By identifying which building footprints intersect with the blackout mask, the analysis isolates the homes likely affected by outages. Converting buildings to centroids allows for precise spatial calculations and comparisons between impacted and non-impacted houses, resulting in an estimated total of approximately 157,970 homes without power.
# Keep buildings that intersect with blackout areas
houses_blackout <- st_intersects(houses, blackout_far)
# Check if each building is in a blackout area
in_blackout <- lengths(houses_blackout) > 0
# Filter buildings in blackout areas
blackout_houses <- houses[in_blackout, ]
# Number of impacted buildings
n_blackout_houses <- nrow(blackout_houses)
# Convert all buildings to centroids
houses_points <- houses %>%
st_centroid()
# Houses impacted by blackouts
houses_impacted <- blackout_houses %>%
st_centroid()
# Houses not impacted by blackouts
houses_not_impacted <- houses_points %>%
filter(!osm_id %in% blackout_houses$osm_id)Let’s create a table to show the number of houses affected by the power outage, making it easier to visualize the impact.
# Table of homes that experienced blackouts
houses_df <- houses_impacted |>
group_by(type) |>
summarise(count = n(),
percent = round((n() / nrow(houses_impacted)) * 100, digits = 2)) |>
st_drop_geometry()
# Add total row
houses_df <- houses_df |>
add_row(type = 'Total', count = sum(houses_df$count))
# Change type NA to character string
houses_df$type[6] <- "NA"
# Display kable table
options(knitr.kable.NA = "")
kbl(houses_df,
table.attr = 'class="custom-table"',
col.names = c(" Types", "Count", "Percentage (%)"),
align = 'c')| Types | Count | Percentage (%) |
|---|---|---|
| apartments | 1136 | 0.72 |
| detached | 353 | 0.22 |
| house | 19760 | 12.51 |
| residential | 1395 | 0.88 |
| static_caravan | 80 | 0.05 |
| NA | 135243 | 85.61 |
| Total | 157967 |
Plot map for the houses that were impacted by the power outage.
# Create the map
tm_basemap('OpenStreetMap') +
tm_shape(houses_not_impacted) +
tm_dots(fill = 'midnightblue', size = 0.02, fill_alpha = 0.8) +
tm_shape(houses_impacted) +
tm_dots(fill = 'firebrick', size = 0.02, fill_alpha = 0.8) +
tm_title('Houses in Houston that lost power',
position = tm_pos_out('center', 'top'),
fontface = 'bold', size = 1.2) +
tm_layout(legend.outside = TRUE,
legend.outside.position = 'right',
legend.title.size = 1,
legend.text.size = 0.8) +
tm_add_legend(labels = c('Impacted houses', 'Not impacted houses'),
fill = c('firebrick', 'midnightblue'),
col = c('firebrick', 'midnightblue'),
alpha = c(0.2, 0.7, 0.5),
title = 'Legend') +
tm_scalebar(position = c('left', 'bottom'),
text.size = 0.6) +
tm_compass(position = c('right', 'bottom'),
type = '4star',
size = 2)The following process links census tract boundaries with median household income data to analyze the socioeconomic characteristics of neighborhoods affected by the blackout. By joining the income dataset to the census tract geometries, the analysis enables comparisons between impacted and non-impacted areas and supports investigation of whether power outages disproportionately affected higher- or lower-income communities.
Let’s verify CRS. Verifying CRS alignment is essential for accurate spatial analysis, as mismatched projections can lead to incorrect geometry intersections, area calculations, or mapping errors.
Warning: Update coordinate reference systems to match!
By aligning the coordinate systems and cropping to the Houston study area, census tracts and blackout affected houses are prepared for spatial comparison. Intersections between the two datasets identify which tracts experienced power outages, and a new column flags these impacted areas. Filtering to only the affected tracts enables a targeted analysis and visualization of how the blackout influenced neighborhoods across the city.
# Transform the CRS and crop to houston bounding box
census_income <- st_transform(census_income, crs = 'epsg:4326')
# Crop census income to Houston
census_income_crop <- census_income %>%
st_crop(houston_bbox)
# Transform and clean blackout houses
blackout_houses_transformed <- blackout_houses %>%
st_transform(crs = 'epsg:4326')
# Identify the blackout areas
census_blackout <- st_intersects(census_income_crop, blackout_houses_transformed)
# Create the column
census_blackout_col <- census_income_crop %>%
mutate(blackout = lengths(census_blackout) > 0)
# Select only the values that are true for plotting
census_blackout_col_filtered <- census_blackout_col[census_blackout_col$blackout, ]Visualize the census tracts in Houston that experienced power outages on a map.
# Map the census tracts that lost power
tm_basemap('CartoDB.VoyagerNoLabels') +
tm_shape(census_blackout_col_filtered) +
tm_polygons(fill = 'blackout',
fill.scale = tm_scale_categorical(
values = c('TRUE' = 'firebrick'),
labels = c('TRUE' = 'Blackout')),
fill.legend = tm_legend(title = 'Blackout status'),
col = 'grey50',
lwd = 0.2) +
tm_title('Census tracts in Houston that lost power',
position = tm_pos_out('center', 'top'),
fontface = 'bold',
size = 1) +
tm_layout(legend.outside = TRUE,
legend.outside.position = 'right',
legend.bg.color = 'white',
legend.bg.alpha = 0.8) +
tm_compass(position = c('right', 'bottom'),
type = 'arrow',
size = 1) +
tm_scalebar(position = c('left', 'top'),
text.size = 0.8)We compare the distributions of median household income between Houston census tracts that experienced blackouts and those that did not during the February 2021 winter storms. By comparing the differences in income, it allows for assessment of whether power outages disproportionately affected higher or lower income neighborhoods, and provides insight into potential socioeconomic disparities in storm impacted areas.
ggplot(census_blackout_col %>% filter(!is.na(B19013e1)), # Remove NA values
aes( x = blackout, y = B19013e1, fill = blackout)) +
geom_boxplot(alpha = 0.8) +
scale_x_discrete(labels = c('TRUE' = 'Blackout',
'FALSE' = 'No blackout')) +
scale_fill_manual(values = c('TRUE' = 'firebrick',
'FALSE' = 'midnightblue'),
labels = c('Blackout', 'No blackout')) +
labs(title = 'Median household income for census tracts',
subtitle = 'Blackout status in Houston census tracts',
x = 'Blackout status',
y = 'Median household income in $',
fill = 'Status') +
scale_y_continuous(labels = scales::dollar_format()) +
theme_minimal() +
theme(plot.title = element_text(face = 'bold', size = 16),
plot.subtitle = element_text(size = 14),
legend.position = 'bottom')This analysis used satellite nighttime light data to identify which parts of Houston lost power during Winter Storm Uri and how these outages related to neighborhood income levels. Surprisingly, the results showed that higher income census tracts appeared slightly more affected, even though past research, including (Lee, Maron, and Mostafavi 2022) found that low income and minority communities experienced more severe and prolonged outages. This difference likely comes from limitations in our dataset. Wealthier neighborhoods often have more outdoor lighting, which can make blackout detection less accurate. Census tract data also hides smaller neighborhood differences, and highway lighting or cloud cover may have caused some areas to be misclassified. In addition, assuming all buildings were residential and using a 200-meter buffer around highways may have excluded areas that actually lost power. Our analysis estimated that about 157,967 residential buildings lost electricity. However, because only two days of satellite data were used, the results may not capture the full timeline of outages. A longer range of dates would make it easier to see where and when blackouts occurred, and provide a clearer picture of storm impacts.
While blackout exposure appeared similar across income levels, the ability to cope with and recover from outages is not equal. Low income communities have fewer resources to stay warm, relocate or rebound after a long outage. This difference between exposure and impact is important for understanding who is most vulnerable during extreme weather events. Future work should use nighttime light data from more days to better track outage duration and recovery. Adding other socioeconomic information, such as age, housing quality, or access to transportation, can improve our understanding of community vulnerability. Combining multiple sources like utility outage records, mobile phone mobility data and satellite imagery would create a more complete picture of who was affected and how. Examining infrastructure characteristics, such as the age of power lines or the placement of substations, may also help explain differences in recovery across neighborhoods. Together, these steps can help identify communities at higher risk and support more effective and equitable preparation for future extreme weather events.
NASA Worldview. (n.d.). NASA Worldview application. NASA Earth Science Data and Information System (EOSDIS). Retrieved October 25, 2025, from https://worldview.earthdata.nasa.gov/
Geofabrik GmbH. (n.d.). Geofabrik OpenStreetMap data extracts. Retrieved October 25, 2025, from https://download.geofabrik.de/
U.S. Census Bureau. (n.d.). American Community Survey (ACS). Retrieved October 25, 2025, from https://www.census.gov/programs-surveys/acs
@online{poudel2025,
author = {Poudel, Aakriti},
title = {Identifying the Impacts of Extreme Weather},
date = {2025-12-03},
url = {https://aakriti-poudel-chhetri.github.io/posts/2025-12-impacts-of-extreme-weather/},
langid = {en}
}